home *** CD-ROM | disk | FTP | other *** search
/ 3D GFX / 3D GFX.iso / amiutils / i_l / irit5 / cagd_lib / cagdoslo.c < prev    next >
C/C++ Source or Header  |  1995-12-30  |  13KB  |  278 lines
  1. /******************************************************************************
  2. * CagdOslo.c - an implementation of the oslo algorithm (1).              *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Aug. 92.                          *
  5. ******************************************************************************/
  6.  
  7. #include <cagd_loc.h>
  8.  
  9. #define KNOT_INFINITY        1e20
  10. #define OSLO_EPSILON        1e-20
  11. #define OSLO_APX_EQ(x, y)    (ABS((x) - (y)) < OSLO_EPSILON)
  12.  
  13. #ifdef DEBUG
  14. static void PrintAlphaMat(BspKnotAlphaCoeffType *A);
  15. #endif /* DEBUG */
  16.  
  17. /*****************************************************************************
  18. * DESCRIPTION:                                                               M
  19. * Computes the values of the alpha coefficients, Ai,k(j) of order k:         M
  20. *                                                                            M
  21. *         n              n    m                   m    n             V
  22. *         _              _    _                   _    _             V
  23. *  C(t) = \ Pi Bi,k(t) = \ Pi \ Ai,k(j) Nj,k(t) = \  ( \ Pi Ai,k(j) ) Nj,k(t)V
  24. *         /              /    /                   /    /             V
  25. *         -              -    -                   -    -             V
  26. *        i=0            i=0  j=0                 j=0  i=0             V
  27. *                                                                            M
  28. * Let T be the original knot vector and let t be the refined one, i.e. T is  M
  29. * a subset of t.                                  M
  30. *   The Ai,k(j) are computed from the following recursive definition:        M
  31. *                                                                            M
  32. *             1, T(i) <= t(i) < T(i+1)                                       V
  33. *           /                                                                V
  34. * Ai,1(j) =                                                                  V
  35. *           \                                                                V
  36. *             0, otherwise.                                                  V
  37. *                                                                            V
  38. *                                                                            V
  39. *           T(j+k-1) - T(i)             T(i+k) - T(j+k-1)                    V
  40. * Ai,k(j) = --------------- Ai,k-1(j) + --------------- Ai+1,k-1(j)          V
  41. *           T(i+k-1) - T(i)             T(i+k) - T(i+1)                      V
  42. *                                                                            M
  43. * LengthKVT + k is the length of KVT and similarly LengthKVt + k is the      M
  44. * length of KVt. In other words, LengthKVT and LengthKVt are the control     M
  45. * points len...                                     M
  46. *                                                                            M
  47. * The output matrix has LengthKVT rows and LengthKVt columns (#cols > #rows) M
  48. * ColIndex/Length hold LengthKVt pairs of first non zero scalar and length ofM
  49. * non zero values in that column, so not all LengthKVT scalars are blended.  M
  50. *                                                                            *
  51. * PARAMETERS:                                                                M
  52. *   k:           Order of geometry.                                          M
  53. *   KVT:         Original knot vector.                                       M
  54. *   LengthKVT:   Length of original knot vector.                             M
  55. *   KVt:         Refined knot vector. Must contain all knots of KVT.         M
  56. *   LengthKVt:   Length of refined knot vector.                              M
  57. *                                                                            *
  58. * RETURN VALUE:                                                              M
  59. *   BspKnotAlphaCoeffType *:   A matrix to multiply the coefficients of the  M
  60. *                              geometry using KVT, in order to get the       M
  61. *                              coefficients under the space defined using    M
  62. *                              KVt that represent the same geometry.         M
  63. *                                                                            *
  64. * KEYWORDS:                                                                  M
  65. *   BspKnotEvalAlphaCoef, alpha matrix, refinement                           M
  66. *****************************************************************************/
  67. BspKnotAlphaCoeffType *BspKnotEvalAlphaCoef(int k,
  68.                         CagdRType *KVT,
  69.                         int LengthKVT,
  70.                         CagdRType *KVt,
  71.                         int LengthKVt)
  72. {
  73.     int i, j, o;
  74.     CagdRType *m, **Rows;
  75.     BspKnotAlphaCoeffType
  76.     *A = (BspKnotAlphaCoeffType *)
  77.                      IritMalloc(sizeof(BspKnotAlphaCoeffType));
  78.  
  79.     A -> Order = k;
  80.     A -> Length = LengthKVT;
  81.     A -> RefLength = LengthKVt;
  82.     A -> Matrix = (CagdRType *) IritMalloc(sizeof(CagdRType) * (LengthKVT + 1)
  83.                                  * LengthKVt);
  84.     Rows = A -> Rows = (CagdRType **) IritMalloc(sizeof(CagdRType *) *
  85.                                   (LengthKVT + 1));
  86.     A -> ColIndex = (int *) IritMalloc(sizeof(int) * LengthKVt);
  87.     A -> ColLength = (int *) IritMalloc(sizeof(int) * LengthKVt);
  88.  
  89.     /* Update the row pointers to point onto the matrix rows. */
  90.     for (i = 0, j = 0; i <= LengthKVT; i++, j += LengthKVt)
  91.     Rows[i] = &A -> Matrix[j];
  92.  
  93.     /* Initialize the matrix with according to order 1: */
  94.     for (m = A -> Matrix, i = (LengthKVT + 1) * LengthKVt; i > 0; i--)
  95.     *m++ = 0.0;
  96.     for (i = 0; i < LengthKVT; i++) {
  97.     CagdRType
  98.         *RowI = Rows[i],
  99.         KVTI = KVT[i],
  100.         KVTI1 = KVT[i + 1];
  101.  
  102.     for (j = 0; j < LengthKVt; j++, RowI++) {
  103.         if (KVTI <= KVt[j] && KVt[j] < KVTI1)
  104.         *RowI = 1.0;
  105.     }
  106.     }
  107.  
  108.     /* Iterate upto order k: */
  109.     for (o = 2; o <= k; o++) {
  110.     for (i = 0; i < LengthKVT; i++) {
  111.         CagdRType
  112.         *RowI = Rows[i],
  113.         *RowI1 = Rows[i + 1],
  114.         KVTI = KVT[i],
  115.         KVTIO = KVT[i + o],
  116.         t1 = KVT[i + o - 1] - KVTI,
  117.         t2 = KVTIO - KVT[i + 1];
  118.  
  119.         /* If ti == 0, the whole term should be ignored. */
  120.         if (t1 < OSLO_EPSILON)
  121.         t1 = KNOT_INFINITY;
  122.         if (t2 < OSLO_EPSILON)
  123.         t2 = KNOT_INFINITY;
  124.  
  125.         for (j = 0; j < LengthKVt; j++, RowI++, RowI1++) {
  126.         int jo = j + o;
  127.  
  128.         *RowI = *RowI * (KVt[jo - 1] - KVTI) / t1 +
  129.             *RowI1 * (KVTIO - KVt[jo - 1]) / t2;
  130.         }
  131.     }
  132.     }
  133.  
  134.     /* Update the row non zero indices. */
  135.     for (i = LengthKVt - 1; i >= 0; i--) {
  136.     for (j = 0; j < LengthKVT && OSLO_APX_EQ(Rows[j][i], 0.0); j++);
  137.     A -> ColIndex[i] = j;
  138.     for (j = LengthKVT - 1; j >= 0 && OSLO_APX_EQ(Rows[j][i], 0.0); j--);
  139.     if ((A -> ColLength[i] = j + 1 - A -> ColIndex[i]) <= 0)
  140.         CAGD_FATAL_ERROR(CAGD_ERR_DEGEN_ALPHA);
  141.     }
  142.  
  143.     return A;
  144. }
  145.  
  146. /*****************************************************************************
  147. * DESCRIPTION:                                                               M
  148. * Frees the BspKnotAlphaCoeffType data structrure.                           M
  149. *                                                                            *
  150. * PARAMETERS:                                                                M
  151. *   A:      Alpha matrix to free.                                            M
  152. *                                                                            *
  153. * RETURN VALUE:                                                              M
  154. *   void                                                                     M
  155. *                                                                            *
  156. * KEYWORDS:                                                                  M
  157. *   BspKnotFreeAlphaCoef, alpha matrix, refinement                           M
  158. *****************************************************************************/
  159. void BspKnotFreeAlphaCoef(BspKnotAlphaCoeffType *A)
  160. {
  161.     IritFree((VoidPtr) A -> Matrix);
  162.     IritFree((VoidPtr) A -> Rows);
  163.     IritFree((VoidPtr) A -> ColIndex);
  164.     IritFree((VoidPtr) A -> ColLength);
  165.     IritFree((VoidPtr) A);
  166. }
  167.  
  168. /*****************************************************************************
  169. * DESCRIPTION:                                                               M
  170. * Same as EvalAlphaCoef but the new knot set NewKV is merged with KVT to     M
  171. * form the new knot vector KVt.                             M
  172. *                                                                            *
  173. * PARAMETERS:                                                                M
  174. *   k:           Order of geometry.                                          M
  175. *   KVT:         Original knot vector.                                       M
  176. *   LengthKVT:   Length of original knot vector.                             M
  177. *   NewKV:       A sequence of new knots to introduce into KVT.              M
  178. *   LengthNewKV: Length of new knot sequence.                                M
  179. *                                                                            *
  180. * RETURN VALUE:                                                              M
  181. *   BspKnotAlphaCoeffType *:   A matrix to multiply the coefficients of the  M
  182. *                              geometry using KVT, in order to get the       M
  183. *                              coefficients under the space defined using    M
  184. *                              KVt that represent the same geometry.         M
  185. *                                                                            *
  186. * KEYWORDS:                                                                  M
  187. *   BspKnotEvalAlphaCoefMerge, alpha matrix, refinement                      M
  188. *****************************************************************************/
  189. BspKnotAlphaCoeffType *BspKnotEvalAlphaCoefMerge(int k,
  190.                          CagdRType *KVT,
  191.                          int LengthKVT,
  192.                          CagdRType *NewKV,
  193.                          int LengthNewKV)
  194. {
  195.     BspKnotAlphaCoeffType *A;
  196.  
  197.     if (NewKV == NULL || LengthNewKV == 0) {
  198.     A = BspKnotEvalAlphaCoef(k, KVT, LengthKVT, KVT, LengthKVT);
  199.     }
  200.     else {
  201.     int LengthKVt;
  202.     CagdRType
  203.         *KVt = BspKnotMergeTwo(KVT, LengthKVT + k,
  204.                    NewKV, LengthNewKV, 0, &LengthKVt);
  205.  
  206.     A = BspKnotEvalAlphaCoef(k, KVT, LengthKVT, KVt, LengthKVt - k);
  207.  
  208.     IritFree(KVt);
  209.     }
  210.  
  211.     return A;
  212. }
  213.  
  214. /*****************************************************************************
  215. * DESCRIPTION:                                                               M
  216. * Prepares a refinement vector for the given knot vector domain with n         M
  217. * inserted knots equally spaced.                         M
  218. *                                                                            *
  219. * PARAMETERS:                                                                M
  220. *   n:       Number of knots to introduce.                                   M
  221. *   Tmin:    Minimum domain to introduce knots.                              M
  222. *   Tmax:    Maximum domain to introduce knots.                              M
  223. *                                                                            *
  224. * RETURN VALUE:                                                              M
  225. *   CagdRType *: A vector of n knots uniformly spaced between TMin and TMax. M
  226. *                                                                            *
  227. * KEYWORDS:                                                                  M
  228. *   BspKnotPrepEquallySpaced, refinement                                     M
  229. *****************************************************************************/
  230. CagdRType *BspKnotPrepEquallySpaced(int n, CagdRType Tmin, CagdRType Tmax)
  231. {
  232.     int i;
  233.     CagdRType dt, t, *RefKV;
  234.  
  235.     if (n <= 0) {
  236.     CAGD_FATAL_ERROR(CAGD_ERR_WRONG_INDEX);
  237.     return NULL;
  238.     }
  239.  
  240.     dt = (Tmax - Tmin) / (n + 1),
  241.     t = Tmin + dt,
  242.     RefKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * n);
  243.  
  244.     for (i = 0; i < n; i++, t += dt) {
  245.     RefKV[i] = t;
  246.     }
  247.     return RefKV;
  248. }
  249.  
  250. #ifdef DEBUG_OSLO
  251. /*****************************************************************************
  252. * DESCRIPTION:                                                               *
  253. * Prints the content of the alpha matrix.                     *
  254. *                                                                            *
  255. * PARAMETERS:                                                                *
  256. *   A:        Alpha matrix to print.                                         *
  257. *                                                                            *
  258. * RETURN VALUE:                                                              *
  259. *   void                                                                     *
  260. *****************************************************************************/
  261. static void PrintAlphaMat(BspKnotAlphaCoeffType *A)
  262. {
  263.     int i, j;
  264.  
  265.     fprintf(stderr, "Order = %d, Length = %d\n", A -> Order, A -> Length);
  266.     for (i = 0; i < A -> Length; i++) {
  267.     for (j = 0; j < A -> RefLength; j++)
  268.         fprintf(stderr, " %9.5lg", A -> Rows[i][j]);
  269.     fprintf(stderr, "\n");
  270.     }
  271.  
  272.     fprintf(stderr, "    ");
  273.     for (j = 0; j < A -> RefLength; j++)
  274.         fprintf(stderr, " %3d %3d |", A -> ColIndex[j], A -> ColLength[j]);
  275.     fprintf(stderr, "\n");
  276. }
  277. #endif /* DEBUG_OSLO */
  278.